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Neutron-star cores contain matter at the highest densities in our Universe. This highly compressed 
matter may undergo a phase transition where nuclear matter melts into deconfined quark matter, 
liberating its constituent quarks and gluons. Quark matter exhibits an approximate conformal 
symmetry, predicting a specific form for its equation of state (EoS), but it is currently unknown 
whether the transition takes place inside at least some physical neutron stars. Here, we quantify this 
likelihood by combining information from astrophysical observations and theoretical calculations. 
Using Bayesian inference, we demonstrate that in the cores of maximally massive stars, the EoS 
is consistent with quark matter. We do this by establishing approximate conformal symmetry 
restoration with high credence at the highest densities probed and demonstrating that the number 
of active degrees of freedom is consistent with deconfined matter. The remaining likelihood is 
observed to correspond to EoSs exhibiting phase-transition-like behavior, treated as arbitrarily rapid 


Strongly interacting matter exhibits deconfined behavior in massive neutron stars 
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crossovers in our framework. 


I. INTRODUCTION 


For macroscopic systems in thermal equilibrium, the 
equation of state (EoS) is a central quantity that re- 
flects not only the basic thermodynamic properties of 
the medium but also its active degrees of freedom 
and thus the underlying phase structure. For neutron 
stars (NSs)—extreme astrophysical objects containing 
the densest matter found in the present-day Universe— 
the EoS is closely related to their measurable macroscopic 
properties due to a link provided by General Relativity 
[1, 2]. This has made NSs a unique laboratory for ul- 
tradense strongly interacting matter, with astrophysical 
observations informing model-agnostic studies of the EoS 
and attempts to determine the phase of matter inside the 
cores of NSs of different masses [3-24]. 

Quantum Chromodynamics (QCD) predicts that at 
very high densities, strongly interacting matter no longer 
resides in a phase where individual nucleons can be iden- 
tified [25-28]. Instead, the active degrees of freedom be- 
come a set of elementary particles, quarks and gluons, 
that form a new phase dubbed quark matter (QM). In 
recent years, ab-initio calculations in nuclear [29-33] and 
particle [34-37] theory have established that the respec- 
tive EoSs of nuclear matter (NM) and QM are qualita- 
tively distinct. While the hadronic EoS is controlled and 
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FIG. 1. Conformalization of neutron-star matter: The 
measure of conformality de = yA? + (A’)?, as a function 
of baryon density. The dark and light red bands correspond 
to 68% and 95% credible intervals (CIs) obtained using a 4- 
segment speed of sound interpolation (ĉa), while the dash- 
dotted lines show the corresponding 68% CIs for a GP re- 
gression. Theoretical limits for de arising from CEFT and 
pQCD are shown as violet bands, while the inferred quantity 
is seen to exhibit a qualitative change in its behavior around 
n ~ 2 —3nsat- The blue, yellow, and green bands show 68% 
credible intervals for the central densities of 1.4Mo,2.0Mo, 
and maximally massive Mrov stars, respectively. The dashed 
horizontal line finally corresponds to our definition of nearly 
conformal matter. The 1.4Mo stars lie firmly above the 
dashed line, while maximally massive stars lie mostly below 
the line. 


characterized by the nucleon mass scale, its QM counter- 
part exhibits near-conformal properties, independent of 
any characteristic mass scales [13, 38]. 

This fundamental difference leaves an imprint on vari- 
ous thermodynamic quantities associated with the equa- 
tion of state, such as the normalized trace anomaly 
A = (e — 3p)/(8¢€), its logarithmic rate of change with 
respect to the energy density A’ = dA/dIne, the poly- 
tropic index y = dInp/dlne, the speed of sound squared 
c? = dp/de, and the pressure normalized by that of a 
system of free quarks, p/pfree. As illustrated in Table 1, 
these quantities indeed take markedly different values in 
(both low- and high-density) NM and in near-conformal 
QM at asymptotically high density. 

At very high baryon densities, where weak-coupling 
methods produce reliable results, strongly interacting 
matter displays near-conformal properties, with confor- 
mality only mildly broken by subdominant loop effects 
and by the small values of the up, down, and strange 
quark masses. As summarized in Table 1 and discussed 
in detail in [13, 38, 39], the characteristics of this sys- 
tem include a small positive normalized trace anomaly 
A, a similarly small but negative A’, a sound speed just 
below its conformal value, c? < 1/3, and a polytropic 
index approaching its conformal limit from above, y = 1. 
Such values are common to many ultrarelativistic sys- 
tems, but in stark contrast with those of the hadronic 
phase, also summarized in Table 1. In the latter case, be 
it at low densities where robust ab-initio results are avail- 
able [3, 29] or at higher densities where one must resort 
to phenomenological models (see, e.g., [13, 40]), the prop- 
erties of hadronic matter are dominated by the nucleon 
mass scale that strongly breaks conformal symmetry. At 
densities exceeding approx. 3n, there is no longer agree- 
ment between different model predictions, but as we dis- 
cuss in detail under Methods, the vast majority of viable 
models remain strongly non-conformal even there. 

It would clearly be very interesting to contrast the 
above expectations with a model-agnostic inference of 
the NS-matter EoS, informed by ab-initio theoretical re- 
sults and astrophysical observations. With few recent 
exceptions [24, 41, 42], existing studies of this kind, how- 
ever, suffer from at least one of two limitations: either 
they fail to take into account high-density information 
from perturbative-QCD (pQCD) calculations, recently 
demonstrated to significantly constrain the NS-matter 
EoS down to realistic core densities [24, 43, 44], or they 
implement observational constraints in the form of hard 
cutoffs, limiting the measurements that can be employed. 

In this work, we remedy the above shortcomings by 
generalizing our earlier analyses [6, 13, 20] to a Bayesian 
framework. This enables us to take advantage of alto- 
gether 12 simultaneous NS mass-radius (MR) measure- 
ments (see the Methods section) and make quantita- 
tive statements about the likelihood of a transition from 
hadronic to QM within stable NSs. Our analysis is per- 
formed using two independent frameworks, the paramet- 
ric speed-of-sound interpolation of [13, 20] and the non- 


parametric Gaussian process (GP) regression of [24, 44], 
which utilizes the high-density constraint from pQCD in 
a considerably more conservative fashion. Both meth- 
ods are used to construct a prior for the EoS, connecting 
a low-density result provided by Chiral Effective Field 
Theory (CEFT) [3, 32] to a high-density limit given by 
pQCD [34, 37]. This prior P(EoS), introduced in detail 
under Methods, is then conditioned using a likelihood 
function incorporating astrophysical measurements, 


data|EoS) P(EoS) 
P(data) i 


P(EoS|data) = A (1) 


with P(data|EoS) = ITI#_,P(data;|EoS) corresponding 
to the uncorrelated individual likelihoods of various NS 
measurements, indexed here by 7. These data, re- 
viewed under Methods, include mass measurements of 
the heaviest pulsars [45-49], the LIGO and Virgo tidal- 
deformability constraints from the NS-NS merger event 
GW170817 [50, 51], and a number of individual mass- 
radius measurements using X-ray observations of pulsat- 
ing, quiescent, and accreting neutron stars by the NICER 
and other collaborations [17, 52-57]. As described under 
Methods, we also vary the number of independent param- 
eters in our parametric interpolation to verify that our 
results are not impacted by an overly restricted prior. 
Additionally, we perform one analysis with a polytropic 
construction to further assess the robustness of our ap- 
proach. The main results from our work take the form 
of posterior distributions for different physical quantities 
evaluated at the centers of NSs of various masses, which 
we compare to theoretical expectations. 


II. RESULTS 
Conformality criterion and theoretical expectations 


The approach of an individual quantity towards its 
conformal value provides a necessary but not sufficient 
condition for the restoration of (approximate) conformal 
symmetry in a physical system. In order to establish the 
conformalization of strongly interacting matter inside NS 
cores, one should therefore simultaneously track several 
quantities and in addition ensure that the system stays 
conformal at higher densities. Here, a particularly useful 
pair turns out to be A and A’. They are related to y and 


2 via 
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implying that when both |A| and |A’| are small, y and 
c? also approach their conformal limits 1 and 1/3. In 
addition, a small value of |A’| naturally ensures that A 
will remain approximately constant at higher densities. 
In order to draw a demarcation line between the non- 
and near-conformal regimes, we combine A and A’ into 


CEFT Dense NM Pert. QM CFTs FOPT 
ce <1 = (0.25, 0.6] <1/3 1/3 0 
A 1/3 [0.05,0.25] — 0, 0.15] 0 1/3-ppr/e 
aN x0 [-0.4,—-0.1] [—0.15,0] 0 1/3- A 
de ~x 1/3  [0.25,0.4] < 0.2 0 >1/(3V2) 
y z 2.5 [1.95, 3.0] [1, 1.7] 1 0 
p/Ptree < 1 (0.25, 0.35] (0.5, 1] —— PPT /Ptree 
TABLE 1. Characteristic features of dense strongly 


interacting matter: Values of a set of six dimensionless 
quantities characterizing the properties of dense strongly in- 
teracting matter in five different limits: “CEFT”, referring to 
predictions for sub-saturation-density nuclear matter [3, 29]; 
“Dense NM” referring to typical nuclear-matter model pre- 
dictions at the highest density (n ~ 3n;) where most models 
still agree with each other, roughly corresponding to the cen- 
ters of typical 1.4Mo pulsars; “Pert. QM” referring to pQCD 
calculations at n = 40nsat with nsat © 0.16fm~? the nuclear 
saturation density [34, 37]; “CFTs” referring to conformal 
field theories in 3+1 dimensions; and “FOPT” referring to a 
system undergoing a first-order phase transition. The indi- 
cated intervals should be considered approximate, while the 
symbol “—” implies the absence of any constraints and ppr in 
the FOPT column refers to the (constant) value of the pres- 
sure at the phase transition. 


a single quantity, adopting the criterion 


de = VA?+(A’)? < 0.2 (3) 
for the identification of near-conformal matter at a given 
density. We justify the value 0.2 as follows. First, as 
noted in Table 1, it represents a natural choice between 
the values this parameter takes in NM and conformal sys- 
tems. Second, at a discontinuous first-order phase tran- 
sition (FOPT), where c? = y = 0 and A’ = 1/3 — A (see 
Eq. (2)), the quantity de can be shown to be bounded 
from below by 1/(3V2) ~ 0.236, so that our criterion pre- 
vents FOPTs from masquerading as conformalized mat- 
ter. We choose to round this value down to 0.2 to pro- 
vide numerical tolerance for the approximate version of 
FOPTs (i.e. arbitrarily rapid crossovers) considered in 
our analysis. Finally, as shown in Table 1, this corre- 
sponds approximately to the largest value obtained for 
this quantity in pQCD. We note that the specific value 
0.2 is a choice, but that our qualitative conclusions are 
insensitive to its small variations. 

While QM is near-conformal, not all matter that is 
near-conformal is QM. Therefore, establishing conformal- 
ization of matter in the cores of NSs does not, in prin- 
ciple, imply the presence of the deconfined phase. One 
quantity not fixed by conformal symmetry is the pres- 
sure normalized by the free (non-interacting) pressure 
P/Pfree. Its value is related to the effective number of 
active degrees of freedom (particle species) Neg in both 
weakly and strongly coupled systems. At weak coupling, 
Dalton’s law states that the pressure is the sum of par- 
tial pressures, which establishes the proportionality of 
Nes = N;N-p/Ptree- This relation can be also derived in 
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FIG. 2. Density dependence of neutron-star-matter 
properties: The normalized trace anomaly A, polytropic 
index y, and speed of sound squared c2 as functions of (left) 
the baryon number density n and (right) the stellar mass M 
normalized by Mrov. The dark and light red bands, dash- 
dotted line and colored bands carry the same meaning as in 
Fig. 1. The CA and GP methods show good agreement for all 
quantities below the maximal density reached in stable NSs, 
and all three quantities display a clear change in behavior 
between the core densities of 1.4Mo and Mrov stars. 


strongly coupled conformal field theories (CFTs) [58, 59] 
as well as in QCD at high temperatures [60, 61]. Because 
this quantity is sensitive to the number of degrees of free- 
dom but insensitive to the interaction strength, we expect 
that QM—be it weakly or strongly coupled—will have to 
exhibit a slowly varying p/Pfree of order one. To this end, 
we may use its behavior to distinguish QM from other 
possible near-conformal phases at NS-core densities. 
From Table 1, we see that in high-density pQCD mat- 
ter p/Pfree is reduced from unity to approximately 0.6 
through perturbative corrections. The normalized pres- 
sure has been extensively studied also in the context 
of high-temperature quark-gluon plasma (QGP), where 
it has played a key role in establishing that deconfined 
matter has been successfully produced in heavy-ion col- 
lisions (see, e.g., [60-62]). Nonperturbative lattice sim- 
ulations have shown that above the pseudo-critical tem- 
perature, its value saturates to a constant around 0.8. 
In other high-temperature quantum field theories in the 
strongly-coupled limit, p/pfree often takes fractional val- 
ues smaller than one [59, 63-67], and for instance in 
N = 4 Super Yang-Mills theory at infinite ’t Hooft cou- 
pling, p/Pfree = 3/4 [59]. In certain lower-dimensional 
CFTs, the value of this quantity is moreover related to 
the central charge of the theory that counts the active 
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properties: Posterior distributions for the polytropic index 
y, speed of sound squared c?, and normalized trace anomaly 
A in the centers of NSs of different masses. The shaded grey 
region corresponds to our conformal criterion de < 0.2, and 
the two-dimensional distributions on the off-diagonal panels 
show 68% and 95% CIs. The one-dimensional distributions 
on the diagonal show the PDFs of the three quantities. 


degrees of freedom [58, 68]. 


Posterior distributions 


The main results of our analysis are displayed in Fig. 1, 
Fig. 2, Fig. 3, and Fig. 4, where we show the behav- 
ior of a number of physical quantities as functions of 
baryon number density, baryon chemical potential, or 
stellar mass. Starting from Fig. 1 and Fig. 2, we show 
credible intervals (CIs) for de and the triplet c2, y, and 
A as functions of the baryon number density n, obtained 
using both a four-segment speed-of-sound interpolation 
(c24) and GP regression. Note that for the c? interpola- 
tion, we display results obtained with a four-segment in- 
terpolation which we have found to be the optimal choice, 
being sufficiently versatile to describe complex EoS be- 
haviors (see Fig. 5 in Methods for results with a varying 
number of segments) and still computationally manage- 
able. 

Concretely, using the a 4 ensemble with the very con- 
servative criterion introduced in Eq. (3), we find the pos- 
terior probability for conformalized matter being present 
in maximally massive stars to be approximately 88%, 
while the corresponding figures for 2Mọ and 1.4Mọo NSs 
are only 11% and 0%, respectively. The same probabil- 
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FIG. 4. Normalized pressure of neutron-star matter: 
The pressure p normalized by the free pressure Pfree as a func- 
tion of baryon chemical potential u, proportional to the ef- 
fective number of degrees of freedom in the system (see main 
text). The dark and light red bands, dash-dotted line, and 
colored bands carry the same meaning as in Fig. 1. This 
quantity is seen to exhibit a plateau within maximally mas- 
sive NSs, reminiscent of an effective restoration of conformal 
symmetry, with a value similar in magnitude to that of weakly 
coupled QM. 


ity for maximally massive stars obtained using the non- 
parametric GP interpolation is 75%, which should be 
considered a robust lower limit for the quantity given 
the very conservative way this method handles the high- 
density constraint. We note that these results are quan- 
titatively stable between different computational choices 
(see Methods), and that they reflect the very conserva- 
tive nature of the new criterion. Had we instead used the 
y < 1.75 criterion introduced in [13], the posterior prob- 
ability of conformalization in maximally massive stars 
would have been 99.8% for the c24 and 97.8% for the GP 
ensemble (see also [41]). Finally, we note that the GP 
method allows a quantitative examination of the impact 
of the pQCD constraint on our results, a topic covered 
more extensively in [24, 69]. As seen in Fig. 6 under 
Methods, this impact is substantial, reflected in the like- 
lihood of QM cores in TOV stars dropping from 75% to 
50% should one not include pQCD input in the analysis. 

In Fig. 3, we next show a detailed investigation of 
the posterior distributions of the quantities displayed in 
Fig. 2. The results are shown both as two-dimensional 
CIs for the joint distributions between pairs of quantities 
and as one-dimensional probability density functions, all 
displayed for both the c24 and GP ensembles. Shown in 
grey is finally the region corresponding to our conformal 
criterion de < 0.2, which is again seen to be considerably 
more restrictive than the criterion y < 1.75 used in pre- 
vious works [13, 41]. For y and c?, we find good quantita- 
tive agreement with previous results, with the inner cores 
of maximally massive stars exhibiting conformalized be- 


havior with high credence. For all three quantities, our 
results as functions of M/Mroy are moreover in good 
agreement with recent works [22, 38, 39, 42, 69, 70], with 
the approach towards conformal behavior as M > Mrov 
clearly visible. Finally, from the y vs. c2 panels, we ob- 
serve that the majority of probability weight for Mrov 
stars that lies outside the grey region resides at small val- 
ues of y and c2. This is indicative of FOPT-like behavior, 
implying that should the cores of Mrovy stars not con- 
tain conformalized matter, then they are destabilized by 
a FOPT, consistent with observations reported in [13]. 


Having established the likely conformalization of mat- 
ter within the cores of maximally massive NSs, we move 
to an analysis of the number of active degrees of freedom. 
In Fig. 4, we show the behavior of the pressure normal- 
ized by its non-interacting limit, p/prree, as a function 
of baryon chemical potential u. Inspecting this result, 
we see a clear flattening of the normalized pressure in 
the vicinity of the cores of maximally massive NSs, with 
a value p/DPfree = 0.40 + 0.03 with 68% credence (corre- 
sponding to Neg = 3.6 + 0.3) that is about 2/3 of the 
central pQCD value in Table 1. That the value of p/Ptree 
within maximally massive NSs is both so close to that 
of weakly coupled QM and appears to evolve very slowly 
gives us confidence in labeling this phase deconfined QM, 
although it should be noted that there exist a number 
of hadronic models, where p/Pfree obtains similar values 
close to the TOV density [40]. As we discuss in Methods, 
in a large majority of the latter models, the quantity is 
already decreasing at this density, which is at odds with 
the behavior we witness. It is, however, interesting to ask 
whether the slightly reduced value of the quantity within 
maximally massive NSs may signal that matter in these 
stellar cores features strongly-coupled characteristics. 


Finally, let us briefly note two interesting issues that go 
somewhat beyond the scope of our current work. First, 
in [71] it was recently suggested that the posterior prob- 
ability of QM cores sharply decreases somewhat below 
Mrov ~% 2.3Mọo, so that if NSs more massive than this 
limit are discovered, QM cores can be ruled out. We have 
failed to reproduce this effect in our results, but instead 
find the probability of QM cores to have only a mild de- 
pendence of Mrov. Second, some explicit FOPTs begin- 
ning at densities below 2ngaz have recently been shown 
to have a noticeable impact on hard limits for the EoS 
[72]. While a generalization of our calculation to include 
explicit discontinuities in the EoS is clearly an interesting 
topic for further work, we again note that the criterion for 
near-conformality we have adopted in this work is very 
conservative. Moreover, the EoSs that [72] identified to 
extend beyond previous hard limits were found to have 
very low posterior probabilities and lie right at the edge of 
the 90% credibility bound for the tidal-deformability con- 
straint. We hence expect that including such EOSs with 
explicit FOPTs within an analysis such at our present 
one would result in only mild shifts in the posteriors. For 
this reason, we trust that the estimates we have obtained 
for the probability of finding conformalized matter inside 


massive NSs are reliable. 


III. DISCUSSION 


In the context of heavy-ion collisions, the properties 
of the equation of state (EoS) have played a key role 
in establishing the production of deconfined hot quark- 
gluon plasma. While there is no discontinuous transi- 
tion in this case, the two phases can still be clearly de- 
marcated based on their qualitatively differing material 
properties. At low temperatures, the EoS exhibits the 
qualitative features of hadronic matter while at high T 
it is characterized by approximate conformal symmetry 
and a number of active degrees of freedom correspond- 
ing to deconfined quarks and gluons. Note, however, 
that approximate conformal symmetry does not imply 
weak coupling: at the temperatures reached in collider 
experiments, weak-coupling methods remain poorly con- 
vergent, and the system exhibits strongly-coupled behav- 
ior. Nevertheless, the transition to QGP is identifiable 
from the restoration of the conformal symmetry, closely 
related to that of chiral symmetry. 

While the high-density matter in NSs in general sig- 
nificantly differs from high-temperature QGP, the pat- 
tern of conformal-symmetry restoration can be similarly 
used as an indicator of phase change inside NSs. In 
this Letter, we have established that the cores of the 
most massive neutron stars are very likely characterized 
by approximate conformal symmetry. We see signs of 
conformalization across multiple microscopic properties 
characterizing the EoS, including the speed of sound c?, 
polytropic index y, normalized trace anomaly A, and its 
logarithmic derivative (with respect to energy density) 
A’. Moreover, we have introduced a new measure for the 
degree of conformality, de = \/A? + (A’)?, the behavior 
of which exhibits a striking transition between 1.4Mo 
and maximally massive NSs, and we have determined 
that the number of active degrees of freedom appears to 
be consistent with theoretical expectations for the prop- 
erties of QM. Concretely, we find that within the cores of 
maximally massive NSs, the ratio p/Pfree is about 2/3 of 
its value in weakly coupled quark matter, leading us to 
the conclusion that the most massive neutron stars very 
likely harbor cores of deconfined QM. An 88% probabil- 
ity (or 75% with the GP method), obtained with a very 
conservative definition of near-conformality, provides a 
quantitative estimate of the remaining uncertainties, of 
which the most important one is related to the possible 
presence of a destabilizing first-order phase transition. 

To conclude, we note that in the context of heavy- 
ion collisions, evidence for the successful production of 
deconfined matter consisted of multiple arguments, of 
which the behavior of the EoS was but one. Similarly, 
it would be desirable to find other observables and ar- 
guments that may either support or contradict the evi- 
dence provided by our analysis of the EoS; for concrete 
suggestions, see e.g. [73-78]. Indeed, just as in the case 


of heavy-ion collisions, to compellingly establish the exis- 
tence of QM in the cores of massive NSs, we need multiple 
lines of evidence all supporting the same picture. 


METHODS 


Our different prior EoS ensembles are constructed by 
joining together three different pieces: a low-density 
CEFT result applied up to approx. 1.2nsat, an interme- 
diate density interpolator, and a high-density part ob- 
tained from pQCD that is used either from u = 2.6 GeV 
onward or robustly translated down to 10nsat using the 
results of [43]. Here, the CEFT part depends of alto- 
gether five parameters in a way specified below, while 
the pQCD pressure only depends on one free parameter, 
the renormalization scale A of the MS scheme. Finally, in 
the intermediate-density region we employ two indepen- 
dent approaches, one built with piecewise-defined poly- 
tropes or the speed-of-sound functions of [13] and the 
other based on a non-parametric Gaussian process re- 
gression introduced in [24]. 

In the case of the parametric interpolation, we sam- 
ple the parameter space with a Markov-Chain-Monte- 
Carlo method using the EMCEE sampler [79]. In practice, 
we use the affine-invariant stretch-move algorithm with 
125 parallel “walkers” to generate Markov-Chain-Monte- 
Carlo chains with Ngnain/T ~ 3 x 10° independent sam- 
ples. Here, 7 ~ 2 x 104 is the autocorrelation time, which 
can sometimes be very long due to the high-dimensional 
parameter space. We have tested the numerical conver- 
gence of each calculation by generating sub-samples of 
the chains (with a jackknife re-sampling) and then veri- 
fying that the resulting credibility regions do not change. 
In addition, for the particular case of the four-segment c2 
interpolation we have verified the convergence by gener- 
ating 50% more samples and noting that all the results re- 
main unchanged. Finally, we smooth the posterior distri- 
butions with kernel density estimation using Silverman’s 
rule to estimate the kernel bandwidth. 

Below, we will in turn cover the construction of our 
low-, high- and intermediate-density EoSs, as well as the 
implementation of the astrophysical constraints. After 
this, we perform a detailed quantitative comparison be- 
tween the three frameworks used at intermediate densi- 
ties (speed-of-sound and polytropic interpolation as well 
as Gaussian Processes), while in the last part of Methods, 
we review the results of an analysis of currently available 
model EoSs for high-density NM. 


A. Low-density EoS 


At low densities, up to slightly above the nuclear sat- 
uration density nsaz, the NS-matter EoS can be reliably 
determined using the well-tested machinery of modern 


nuclear theory, including the systematic effective theory 
framework of CEFT [29-33]. Specifically, we take the 
crust EoS from [80], while starting from the crust-core 
transition around 0.5nsaz, we use the results of [32] for 
pure neutron matter (PNM) based on the N3LO ab-initio 
calculations using the A = 500 MeV potentials from [31], 
reported up to 0.195 fm~3, or slightly above 1.2ngat.. 
The upper limit for the pressure is chosen so that all 
(sampled) PNM EoSs have non-negative speeds of sound 
squared throughout the low-density interval considered, 
and we similarly ensure the same condition after beta 
equilibration. 

To transform the PNM results into a beta-stable form, 
we utilize the Skyrme-motivated ansatz introduced in [3], 
which expresses the sum of the interaction and kinetic 
energies in the form 


Tsum 5 
— [2(a — 2az)z(1 — x) + az] n (4) 
+ [2(n — 2nz)æ(1 — x) +n] n". 


alma) 3 [45/9 4 (1 —2)%9] (any 


Here, x = n,/n denotes the proton fraction, i = n/Nsat, 
and a, az, 7, and ny are parameters that can be deter- 
mined from the saturation properties of symmetric nu- 
clear matter (see [3] for discussion and numerical val- 
ues). The normalization constant Tgsnm on the other 
hand stands for the Fermi energy of symmetric nuclear 
matter (SNM, zx = 1/2) at n = neat, 


1 37? Neat 2/3 
T; =e 5 
SNM ~ Smg ( 2 ) (5) 


where mp is the (mean) baryon mass, taken equal to that 
of the neutron my œ% 939.565 MeV. Finally, we add to 
Eq. (4) an additional phenomenological term, 


where we have introduced two further phenomenological 
parameters ño and CL. Their presence can be seen to lead 
to a considerable improvement in the model’s ability to 
describe the CEFT data for PNM. 

Following the original study [3], we next require that 
the energy and pressure satisfy 


1. €(Nsat, 1/2) = —16 MeV, 
2. p(Msat, 1/2) = 0, 


where € = €; + €2 and p = n70e€/On, and which aids us 
in the selection of the model parameters œ and 7. At 
the same time, we do not fix the parameter [ using a 
single value of the incompressibility parameter for SNM 
at saturation 


2 Oe 


Ks = Msat a3 


(Nsat» 1/2) (7) 


as done in [3], but rather consider I a free parameter fol- 
lowing [19]. Moreover, with the mass difference between 
the proton and neutron approximated to zero, we note 
that the proton fraction x = n,/n can be solved from the 
equation 


OO) i 1(F,2) =0, (8) 
t 
where ue = V3722n is the electron chemical potential in 
the ultra-relativistic limit [3]. 

In practice, we draw a large sample of the PNM EoSs 
of [32] using a prior distribution for the five PNM model 
parameters: az, nL, [, Cz, and ño. This distribution 
is then further approximated using a Gaussian mixture 
model to get a simpler and smoother multidimensional 
distribution, which is used to represent the CEFT result 
in building the NS-matter EoSs. 


B. High-density EoS 


Due to the asymptotic freedom of QCD, at very high 
densities the behavior of cold (T = 0) strongly interacting 
matter can be approached from the deconfined side us- 
ing the machinery of modern perturbative thermal field 
theory. Here, we rely on state-of-the-art perturbative- 
QCD (pQCD) calculations for three-flavor quark matter, 
presented in detail in [34, 37]. 

In the limit of beta equilibrium and strictly vanishing 
temperature, the pQCD pressure depends on only two 
parameters: the baryon chemical potential u and the MS 
renormalization scale A. The latter is often scaled by the 
central value of two times the quark chemical potential 
2/3 to create a dimensionless parameter X = 3A/(2y1), 
which is typically varied from 1/2 to 2. For the prior 
probability density on X, we use 


1 
[In(Xmax/Xmin)| X 


f(X) = (9) 
with Xmin = 1/2 and Xmax = 2. The pQCD pres- 
sure is then used for baryon chemical potentials u > 
UpQocp = 2.6 GeV, corresponding to baryon densities 
NpQCD 2 A0Nsat- 


C. Intermediate densities 


As noted in the main text, at intermediate densities 
we use two independent and complementary methods for 
constructing the full EoSs, largely following earlier works 
[13, 20, 24]. As none of the methods we use here is new, 
we leave a more extensive technical introduction to the 
original references. 

In the case of parametric interpolation, we use both 
a piecewise-polytropic or piecewise-linear-c? EoSs with a 
varying number N of intermediate segments between the 
CEFT and pQCD regimes. The resulting EoSs contain 


2N +1 or 2N free parameters in the polytropic and c? 
cases, respectively. 

Finally, for the nonparametric GP calculations, we ex- 
actly follow the prior construction of [24]. 


D. Astrophysical measurements and constraints 


In addition to the theoretical constraints discussed 
above, our EoS models have been further conditioned us- 
ing astrophysical observations. We list and discuss these 
measurements here. 

We first discuss the pulsar mass measurements. We 
demand that an EoS be compatible with the latest mass 
measurements of massive pulsars. So far, two NS sys- 
tems have been observed likely containing a two-solar- 
mass NS, PSR J0348+0432 (M = 2.01 + 0.04 Mo, 
68% confidence interval) [46] and PSR J0740+6620 
(M = 2.08 + 0.07 Mo, 68% credible interval) [49]. We 
model these mass measurements as normally distributed, 
Le. Moa32/Mo ez N (2.01, 0.04?) and Mes20/Mo ma 
N (2.08, 0.077), respectively. 

We then require that the TOV mass derived from a 
target EoS has to be greater than the observed masses 
of these two stars; otherwise, the EoS in question is dis- 
carded. 

We now discuss the tidal deformability measurements. 
We use GW data from the first binary NS merger 
event GW170817 collected by the LIGO/Virgo collab- 
oration [50]. In particular, we use the marginalized 2d 
posterior distribution for the tidal deformabilities of the 
two merger components, given in Fig. 1 of [51]. The indi- 
vidual tidal deformabilities for our EoSs are on the other 
hand calculated using formulae given in Ref. [81] assum- 
ing that both compact objects are NSs. Note that the 
LIGO/Virgo collaboration has also detected another pos- 
sible NS-NS merger event GW190425 [82]. It is, however, 
unfortunately too faint to be of use to constrain the NS- 
matter EoS and is therefore omitted. Moreover, the col- 
laboration has also detected three possible BH-NS merger 
events: GW190814 [83], GW200105, and GW200115 [84]. 
It is, however, very likely that GW190814 is just another 
BH-BH event. 

In our analysis, we use two prior variables, the chirp 
mass of the binary 


(mima) 

M= rema E (10) 
and the mass ratio q = m2/m,, where m, > mo, 
instead of the individual masses mı and mə of the 
binary components. The LIGO/Virgo collaboration 
has determined the binary chirp mass to be M = 
1.186 + 0.001 Mo (90% symmetric credibility lim- 
its) [85], which we take to be normally distributed, i.e. 
M/Mọo ~ N(1.186,3.7 x 1077). The distribution of the 
mass ratio q is on the other hand chosen to be uniform 
between 0.4 and 1, i.e. q ~ U(0.4,1). Finally, we discard 
all parameter candidates that do not satisfy 


1. m1, M2 < Mrov and 


We now turn to the mass-radius measurements from 
X-ray observations. The generated EoSs are also con- 
strained by 12 astrophysical NS mass-radius measure- 
ments that are introduced in Table 2 together with the 
chosen mass priors. The corresponding radius of an indi- 
vidual star can be determined using the TOV equations 
when the EoS is generated and a candidate mass is drawn 
from a uniform prior distribution, m; ~ U(Mmin, Mmax); 
where Mmin and Mmax are the minimum and maximum 
masses, respectively. The limits of each mass prior are 
selected such that they encapsulate the complete range 
of the mass-radius measurement, with a maximal range 
of m;/Mo € (0.5, 2.8]. When possible, tighter limits are 
selected because they lead to more efficient sampling of 
the distribution We also ignore the additional factor to 
the prior, ~ (mroy — Mmin)~1, correcting for the mass 
selection bias, as the factor is to a good approximation 
constant for our choice of Mmin ¥ 0.5Mo [10, 15]. 

Proceeding to the individual measurements considered, 
we first use the 2d mass-radius probability distributions 
of the pulsars PSR J0030+0451 and PSR J0740+6620 as 
reported by the NICER mission. The PSR J0030+0451 
measurement is based on publicly available data from 
[56] employing their ST+PST model, while the PSR 
J0740+6620 measurement is similarly based on publicly 
available data from [17]. From the latter reference, we 
employ their combined NICER+XMM constraints to- 
gether with the inflated cross-instrument calibration er- 
ror. The differences between this choice and other pos- 
sible models would result in small radius differences of 
AR < 0.1km. 

In addition to the above, we also incorporate 2d mass- 
radius measurements from three different X-ray bursters. 
The most accurate constraint corresponds to the neutron 
star in the binary system 4U J1702—429 and is derived 
using direct atmosphere model fits to time-evolving en- 
ergy spectra [54]. Here, we use their model D measure- 
ment corresponding to a fit with a free mass, radius, and 
composition. The other two neutron star mass-radius 
measurements, corresponding to the binary systems 4U 
1724—307 and SAX J1810.8—260 [86], are derived us- 
ing a cooling tail method analysis [87]. We do not con- 
sider other possible bursting sources because of the large 
uncertainty related to the unconstrained accretion disk 
state [88]. We also note that no rotational effects are 
modeled in [89]; these can introduce an additional uncer- 
tainty of order 0.5 km into the radius of a rapidly rotating 
source. 

Lastly, we use mass-radius measurements of seven neu- 
tron stars in quiescent low-mass X-ray binary systems 
[52, 53]. We only use systems with reliable distance mea- 
surements, which enables breaking the degeneracy be- 
tween the observed flux and emitted luminosity (i.e., the 
size of the emitting region and the source distance). We 
use models that assume the surface to be uniformly emit- 


System Mass prior [Mo] Model Ref. 
NICER pulsars 
PSR J0030+0451 (1.0, 2.5) ST+PST [55, 56] 


PSR J0740+6620 N(2.08,0.077) N+XMM+cal. [17, 49, 57] 


qLMXB systems 


M13 U(0.8, 2.4) H 52 
M28 U(0.5, 2.8) He 53 
M30 U(0.5, 2.5) H 53 
w Cen U(0.5, 2.5) H 53 
NGC 6304 Uu(0.5, 2.7) He 53 
NGC 6397 U(0.5, 2.0) He 53 
47 Tuc X7 u (0.5, 2.7) H 53 
X-ray bursters 
4U 1702—429 U(1.0, 2.5) D 54 
4U 1724—307 U(0.8, 2.5) SolA001 86 
SAX J1810.8—260  U(0.8, 2.5) SolA001 86 


TABLE 2. Neutron-star measurements: A summary of 
the NS X-ray M-R measurements considered in this work. 


ting, indicating the absence of hot/cold spots. The atmo- 
sphere composition is selected manually (from H vs. He) 
using a composition that gives a radius of R ~ 12km; 
varying the composition causes a large (~ 2x) change in 
R, so the selection of a realistic composition is relatively 
unambiguous for all of the used sources. 

Finally, we discuss our implementation of multimes- 
senger constraints arising from GW170817. The timing 
and spectral properties of the short gamma-ray burst 
(sGRB) and kilonova from the GW170817 event are also 
used to set additional constraints for the EoS. These 
constraints are based on the astrophysical jet genera- 
tion and launching models that have demonstrated that 
the jet launching from a black hole merger remnant is 
strongly favored for sGRB events. Additionally, the 
prolonged existence of a supramassive/hypermassive NS 
(SMNS /HMNS, respectively) remnant is disfavored be- 
cause of the observed moderate kinetic energy of the kilo- 
nova/sGRB: an SMNS/HMNS inside a thick disk would 
release substantial rotational energy into the surround- 
ing kilonova, strongly modifying its appearance. No 
such modifications, such as a prolonged prompt emis- 
sion phase or strong blue-shifted winds, were detected 
from GW170817. This implies that the merger remnant 
either (i) collapsed immediately into a black hole or (ii) 
formed an SMNS or HMNS remnant that then shortly 
after collapsed into a black hole (see, e.g. [20] and refer- 
ences therein). 

The short gamma-ray burst from GW170817 was de- 
tected with a lag of < 2s after the NS merger, mea- 
sured as the difference between the peak of the gravita- 
tional wave signal and the arrival of the first y-rays. We 
therefore require that the remnant must (at least) form 
a supramassive NS and therefore has a total (baryon) 
mass exceeding the (baryon) TOV mass, mp1 + Me, 2 > 
mM rov. We assume here zero ejecta for the merger pro- 
cess, following previous works (see [20] for references and 


discussion). Additionally, we note that the short lag of 
< 2s favors the HMNS scenario, given that SMNSs are 
expected to be more long-lived with typical lifetimes of 
2 10s. Enforcing the HMNS scenario would lead to an 
even more stringent constraint, Mı + M2 > mgmng, but 
following [20], here we only consider the more relaxed 
constraints following from SMNS formation. 


E. Comparison of the constructed EoSs 


Our fiducial EoS consists of a piecewise-linear c? 


terpolant with N = 4 intermediate segments as a 
function of the baryon chemical potential u, for which 
we use the shorthand notation c24. We have bench- 
marked the robustness of the corresponding results to 
other interpolation methods, including a varying num- 
ber of interpolation segments, with results of this com- 
parison displayed in Fig. 5. Comparing the parame- 
ter regions that the models can sample for de, c?, Y, 
p, and the NS mass-radius relation, we find that four 
segments provides the minimum viable interpolation ac- 
curacy. The biggest shortcoming of the computation- 
ally most-economic three-segment model (i.e., c23) is the 
lack of EoSs that mimic phase transitions at densities of 
n & 5 — 10nsat, whereas only negligible deviations are 
found between the four-segment and five-segment mod- 
els. We note in passing that the polytropic interpolation 
(pn) gives very similar results when the number of seg- 
ments JN is at least four. 


in- 


A similar comparison of the effects of various astro- 
physical observations on our results is shown in Fig. 7 
and Fig. 6, where we separately consider the GP and 
Ga methods, respectively. In particular, we compare 
the resulting Bayesian posterior distributions for calcu- 
lations with no astrophysical measurements (denoted as 
“no obs.”); only mass-measurements from radio pulsars 
(“r”); pulsar mass measurements, GW deformabilities, 
and the SMNS hypothesis for GW170817 (“r + GW”); 
and finally pulsar mass measurements, GW deformabil- 
ities, the SMNS hypothesis for GW170817, and X-ray 
measurements (“r+GW-+X-ray”). Our fiducial calcula- 
tions are based on the measurement set “r+GW-+X-ray” , 
which includes all measurements. We note that espe- 
cially de does not strongly depend on the inclusion of X- 
ray measurements, and that there are some interesting 
quantitative differences between the responses of the ce 4 
and GP ensembles to the astrophysical measurements. 
This may be partially related to the fact that in Fig. 6, 
the high-density constraint from pQCD is only included 
in the fifth additional column, owing to the fact that in 
this approach, it is possible to turn this effect on and 
off at will. Once all observations and the pQCD input 
have been included, our final results for different physical 
quantities obtained using the two methods largely agree. 


F. Properties of nuclear matter models at high 
density 


In this final part of the Methods section, we provide 
brief justification for our claim that viable models of 
dense nuclear matter behave differently from our results, 
reflected in part in the values appearing in the “Dense 
NM” column of Table 1. To study this question in a sys- 
tematic fashion, we have analyzed a large set of publicly 
available nuclear-matter models, determining the charac- 
teristic ranges of various thermodynamic quantities they 
predict. The result of this analysis is given in Fig. 8, 
which display results for all hadronic T = 0 EoSs ap- 
pearing in the CompOSE database [90]. (See also a more 
detailed analysis in the framework of relativistic mean- 
field models [91] with qualitatively similar results for c2, 
y and A.) Note that here we do not discard EoSs that 
fail to satisfy given hard observational limits (unlike in 
the discussion of [13]), but rather impose a color coding 
for the curves corresponding to their relative likelihoods 
in comparison with the most likely EoS in our GP en- 
semble. We also emphasize that the EoSs from the Com- 
pOSE database we analyze are not required to conform 
to our low-density CEFT prior; the latter appears only in 
the priors of our EoS ensembles and not in the likelihood 
function that we use here. 

Our first observation from Fig. 8 (panel a) is that the 
different NM models appear to be in qualitative agree- 
ment with each other for densities n S 3ns, where they 
all display strongly non-conformal behavior: a large or 
rapidly changing trace anomaly, a large value of the poly- 
tropic index y, a rapidly rising value of the speed of sound 
c?, and most indicatively, a large value of the conformal 
distance de. Above this density, the agreement is quickly 
lost (visible in particular in de), which is why the ranges 
we report for different quantities in Table 1 are chosen 
to correspond to n = 3ns. Proceeding towards TOV 
densities, denoted by dots at the end of the curves in 
the four plots, some model predictions even end up at 
near-conformal values, but none of these models carries 
a non-negligible likelihood in our analysis. Discarding 
them, the lower end of the range we give for de in Ta- 
ble 1 is observed to stay valid at all densities, so that 
sizable de values can be considered a characteristic pre- 
diction of NM models. Finally, in panel b of the same 
figure we display the behavior of the normalized pressure 
in the same NM models, noting that in most cases the 
values of the quantity are undergoing a slow decline at 
the TOV point. This appears to be at odds with the 
general trend of the interpolated EoS band. 


DATA AVAILABILITY 


Thinned versions of the EoS ensembles generated for 
this study (interpolations and GP) have been deposited 
to the Zenodo data repository and are available sepa- 
rately for the interpolations [92] and the GP method [93]. 


CODE AVAILABILITY 


All codes written for the generation and analysis of 
the EoS ensembles are available from the authors upon 
request. 
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FIG. 5. Comparison of different interpolators: A comparison of the parametric (interpolated) EoS model results, illus- 
trating their robustness under varying the number of parameters as well as the form of the interpolants. Shown are always 
posterior CIs from the full analysis with all astrophysical constraints. The first column shows the polytropic interpolation with 
four segments (p4). The following columns show the c2 interpolation with an increasing number of linear segments from 3 
(c23) to 4 (c24) to 5 (c25), as well as the Gaussian processes parameterization (GP). For the 5-segment c? ș interpolation, the 
posterior sampling has not yet formally converged and, hence the results are only indicative. The rows show the evolution of 
the conformal measure de, the speed of sound squared c2, polytropic index y, mass-radius M-R relation, and pressure-energy- 
density p-e. The colorbars are tailored specifically to best highlight the variations in the probability density. 


GP (no obs.) 


GP (r + GW) 


0.65 
20.44 
vV 


0.24 


1 


GP (r + GW + X-ray) GP (r + GW + X-ray + QCD) 


i] 
l] 
i 
i 
1 
1 


0° n [Msat] 


o” n [Nsat] 


f 
10! 10° n [nsa] 


10! 10° 


n [Nsat] 


n [Msat] 


10! 10° 


n [Nsat] 


10! 10° n [nsa] 


10! 10° 


n [Nsat] 


10! 10° n [nsa] 


rrit 
10! 10° 


n [Msat] 


rrr 
10! 10° 


n [Nsat] 


10° n [Msat] 


riririlyprsi divi tiiiitinn gg 


12 13 
R [km] 


p [MeV fm? 


e [MeV fm™?] 


FIG. 6. 


T, 
10° 


e [MeV fm7>] 


10° 


e [MeV fm™?] 


rane 
10° 


e [MeV fm7?] 


e [MeV fm™?] 


12 


Effect of different measurements — Gaussian Processes: Model comparison showing the effect of different 


measurements on the GP parameterization. The first column shows the CIs of various quantities with an EoS that does not 
have any astrophysical measurements nor high-density QCD input (prior). The second column shows the Bayesian posterior 
densities for an EoS conditioned with the ~ 2Mo radio pulsar mass measurements. The third column shows the calculations 
with pulsar masses, tidal deformabilities from GW170817, as well as the assumption that the remnant in GW170817 collapsed 
into a black hole. The fourth column shows the calculation with pulsar masses, tidal deformabilities, the assumption that 
GW170817 collapsed to a BH, and X-ray measurements. The fifth column shows the calculations with all the astrophysical 
input from the fourth column, along with the addition of the high-density pQCD input. The gray vertical bands correspond 
to the density interval found at the centers of TOV-mass stars (with 68% confidence level). Other quantities and symbols as 


in Fig. 5. 
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FIG. 7. Effect of different measurements — interpolation: Model comparison showing the effect of different measure- 
ments on the ea parameterization. The first column shows the CIs of various quantities with an EoS parameterization that 
does not have any astrophysical measurements (prior). The second column shows the Bayesian posterior densities for EoS 
conditioned with the ~ 2Mo radio pulsar mass measurements. The third column shows the calculations with pulsar masses 
and tidal deformabilities from GW170817, as well as the assumption that a BH was formed in GW170817. The fourth column 
shows the calculation with pulsar masses, tidal deformabilities, the assumption that a BH was formed in GW170817, and X-ray 
measurements. Other quantities and symbols are as in Fig. 5. 
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FIG. 8. Comparison to nuclear-matter models: Panel a: The physical quantities displayed in Fig. 1 and Fig. 2, but 
this time including as the thin lines all purely hadronic (T = 0, 8-equilibrium) models available on the CompOSE database 
[90]. The coloring of the individual EoSs corresponds to the relative likelihoods of the models in comparison with the most 
likely EoS from the GP ensemble. This likelihood depends only on the astrophysical and high-density pQCD inputs. The 
pQCD constraint is imposed at the TOV points where these curves end, and is implemented by testing whether these points 
can be connected to the pQCD limit at 40n, in a causal and thermodynamically consistent fashion. The solid red vertical 
lines here represent the ranges displayed in Table 1, and the dark and light red regions in the background correspond to our 
Gea interpolation results. Panel b: The normalized pressure of Fig. 4, but including now also the purely hadronic (T = 0, 
8-equilibrium) model EoSs displayed in panel a. The curves are again set to end at the respective TOV points. 
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